/* Copyright (C) 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/***********************************************************
 * lammps.c -- functions for interfacing Brahms and LAMMPS * 
 ***********************************************************/

#include "dataStructs.h"

extern const VecR region, *siteCoords, *siteOrientVec; 
extern const Site *site;
extern const real rCutLipLip;
extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps, massAmine, scaleRefSpringLength, refAngle_orient;
extern const real springRigidity, angleRigidity, orientRigidity, refAngle_CholPhosGly, refAngle_PhosGlyEst, refAngle_cisUnsat;
extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters, excludeBondedNebrs;
extern int zReplicas;

// local variables:
static const int nAtomsPerLipid = 15; // assuming 15-site lipid models
static const int nHarmonicBondsPerLipid = 14; // assuming 15-site lipid models

static const int nBondTypes = 5; // types defined below: 
static const int X_PHOS = 1;
//static const int CHOL_PHOS = 1;
//static const int AMINE_PHOS = 7;
static const int PHOS_GLY = 2;
static const int GLY_EST =3;
static const int EST_CHAIN = 4;
static const int CHAIN_CHAIN = 5;

static const int nCosineAnglesPerLipid = 13; // assuming 15-site lipid models: 13
static const int nDipoleAnglesPerLipid = 3; // assuming 3 dipoles to restrain
static const int nAngleTypes = 5; // 1_chol(amine)-phos-gly, 2_phos-gly-est, 3_PI_sat, 4_cisUnsat, 5_dipoles
static const int _phos_ = 1;
static const int _gly_ = 2;
static const int _piSat_ = 3;
static const int _cisUnsat_ = 4;
static const int DIPOLE = 5;

/*******************************************************************************************************
 * WriteLammpsData -- writes "data." file to be read into the LAMMPS program through read_data command *
 *******************************************************************************************************/
void WriteLammpsData( FILE *fp )
{
  int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3, zRep;
  int idGly, idEstA, idEstB, idEst, idChain1, idChainI, idChainJ;
  int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS topology
  int nBonds = nHarmonicBondsPerLipid*nLipids;
  int nAngles = (nCosineAnglesPerLipid+nDipoleAnglesPerLipid)*nLipids;
  int nAtomTypes;

  if ( nLipids == nDOPCsDSPCs || nLipids == nDOPEs ) nAtomTypes = nTypes-1; // bilayer containing only one lipid species
  else nAtomTypes = nTypes; // mixed bilayer

  if ( zReplicas==0 ) zReplicas = 1; // set zReplicas to default value in case it wasn't set in 'brahms.md'

  nAtoms *= zReplicas;
  nBonds *= zReplicas;
  nAngles *= zReplicas;

  fprintf(fp, "LAMMPS 'data.' description - created with Brahms (www.personal.soton.ac.uk/orsi/brahms)\n");
  fprintf(fp, "\n");

  fprintf(fp, "%9d atoms\n", nAtoms);
  fprintf(fp, "%9d bonds\n", nBonds);
  fprintf(fp, "%9d angles\n", nAngles);
  fprintf(fp, "\n");
  fprintf(fp, "%9d atom types\n", nAtomTypes);
  fprintf(fp, "%9d bond types\n", nBondTypes);
  fprintf(fp, "%9d angle types\n", nAngleTypes);
  fprintf(fp, "\n");

  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.x*LENGTH_Angstrom/2, region.x*LENGTH_Angstrom/2, "xlo xhi");
  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.y*LENGTH_Angstrom/2, region.y*LENGTH_Angstrom/2, "ylo yhi");
  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.z*LENGTH_Angstrom/2, region.z*(zReplicas-0.5)*LENGTH_Angstrom, "zlo zhi");
  fprintf(fp, "\n");
  
  fprintf(fp, "Atoms\n\n");
  atomCount = 0;
  for ( zRep = 0; zRep < zReplicas; zRep++ ) {
    DO_SITE {
      atomCount++;
      if ( site[n].type==AMINE_TYPE )
	fprintf(fp, "%7d%5d", atomCount, 2 ); 
      else fprintf(fp, "%7d%5d", atomCount, site[n].type+1 ); // write atom-ID and atom-type
      fprintf(fp, "%11.4f%11.4f%11.4f", siteCoords[n].x*LENGTH_Angstrom, siteCoords[n].y*LENGTH_Angstrom, (siteCoords[n].z+zRep*region.z)*LENGTH_Angstrom); 
      if ( site[n].inLipid == 0 ) // water 
	fprintf(fp, "%7d", site[n].inLipid); // write molecule-ID (this is 0 for water sites, 1 for lipid-1, 2 for lipid-2, etc.)
      else fprintf(fp, "%7d", site[n].inLipid+(zRep*nLipids) ); // write molecule-ID (this is 0 for water sites, 1 for lipid-1, 2 for lipid-2, etc.)
      fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this is 0.00 for point-dipole sites):
      fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole components (non-zero only for water, glycerol and ester sites)
	      siteOrientVec[ n ].x * dipole[site[n].type] * DIPOLE_D * D_IN_eA, 
	      siteOrientVec[ n ].y * dipole[site[n].type] * DIPOLE_D * D_IN_eA, 
	      siteOrientVec[ n ].z * dipole[site[n].type] * DIPOLE_D * D_IN_eA );
      // write attributes for 'sphere' style - diameter & atom density:
      fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom, mass[site[n].type] / (4*PI/3 * Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
      fprintf(fp, "\n"); 
    }
  }
  if (atomCount!=nAtoms) {printf("Error: atomCount = %d != %d = nAtoms.\n", atomCount, nAtoms); exit(0);}
  fprintf(fp, "\n"); 
  
  fprintf(fp, "Bonds\n\n");
  bondCount = 0;
  for ( zRep = 0; zRep < zReplicas; zRep++ ) {
    DO_SITE {
      id = (n + 1) + zRep * ( nAtoms / zReplicas ); // conversion to LAMMPS atom-ID
      switch ( site[ n ].lipidUnit ) {
      case WATER:     /* break - nonlipid site */ 
      case TAIL_A5: // break - this case is considered in TAIL_A4
      case TAIL_B5: // break - this case is considered in TAIL_B4
	break;
      case CHOLINE_OR_AMINE: /* choline***phosphate (DOPC) or amine***phosphate (DOPE) */
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, X_PHOS, id, id+1 );
	break;
      case PHOSPHATE: // phosphate***glycerol 
	atom1 = id; // phos
	atom2 = atom1 + 1; // gly
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, atom1, atom2 ); 
	break;
      case GLYCEROL: 
	idGly = id; // gly
	idEstA = idGly + 1; 
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstA ); // GLYCEROL***ESTER_A:
	idEstB = idGly + 7; // because there are 7 sites between glycerol and esterB in the topological definition of a lipid molecule
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstB ); // GLYCEROL***ESTER_B:
	break;
      case ESTER_A: // ESTER_A***TAIL_A1
      case ESTER_B: // ESTER_B***TAIL_B1
	idEst = id;
	idChain1 = idEst + 1; 
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, EST_CHAIN, idEst, idChain1 ); 
	break;
      case TAIL_A1: // TAIL_A1***TAIL_A2
      case TAIL_A2: // TAIL_A2***TAIL_A3
      case TAIL_A3: // TAIL_A3***TAIL_A4
      case TAIL_A4: // TAIL_A4***TAIL_A5
	idChainI = id;
	idChainJ = idChainI + 1;
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI, idChainJ ); 
	break;
      case TAIL_B1: // TAIL_B1***TAIL_B2
      case TAIL_B2: // TAIL_B2***TAIL_B3
      case TAIL_B3: // TAIL_B3***TAIL_B4
      case TAIL_B4: // TAIL_B4***TAIL_B5
	idChainI = id;
	idChainJ = idChainI + 1;
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI, idChainJ ); 
	break;
      default:
	printf("Error: default case reached in WriteLammpsData. Exiting...\n");
	exit(0);
      }    
    }
  }
  if (bondCount!=nBonds) {printf("Error: bondCount=%d!=%d=nBonds.\n",bondCount,nBonds); exit(0);}
  fprintf(fp, "\n");

  fprintf(fp, "Angles\n\n");
  angleCount = 0;
  for ( zRep = 0; zRep < zReplicas; zRep++ ) {
    DO_SITE {
      id = n +  zRep * ( nAtoms / zReplicas ); // conversion to LAMMPS atom-ID
      switch ( site[ n ].lipidUnit ) { /* identify the site "n" corresponding to the angle "(n-1)--(n)--(n+1)" */
      case WATER: 
      case CHOLINE_OR_AMINE:
      case TAIL_A5: 
      case TAIL_B5: 
	break; /* there is no angle potential centered on these sites */
      case PHOSPHATE: /* choline_or_amine(1)-phosphate(2)-glycerol(3) */
	atom1 = id; // choline or amine
	atom2 = id+1; // phosphate
	atom3 = id+2; // glycerol
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, atom1, atom2, atom3 ); 
	break;
      case GLYCEROL:
	atom1 = id; // phosphate
	atom2 = id+1; // glycerol
	atom3 = id+2; // ester_a
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3 ); // phos-gly-estA
	atom3 = id+8; // ester_b
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3 ); // phos-gly-estB
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, atom2, atom1, atom1 ); // gly-phos angle_dipole
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, id+2, id+3, id+3 ); // estA-tailA1 angle_dipole
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, id+8, id+9, id+9 ); // estB-tailB1 angle_dipole
	break;
      case ESTER_A: /* glycerol(3)-ester_A(4)-hydrocarbon_a1(5): */ 
      case TAIL_A1: /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */  
      case TAIL_B1: /* ester_B(4)-hydrocarbon_b1(5)-hydrocarbon_b2(6) */
      case TAIL_A3:   /* hydrocarbon_a2(6)-hydrocarbon_a3(7)-hydrocarbon_a4(8) */
      case TAIL_A4:   /* hydrocarbon_a3(7)-hydrocarbon_a4(8)-hydrocarbon_a5(9) */
      case TAIL_B3:   /* hydrocarbon_b2(6)-hydrocarbon_b3(7)-hydrocarbon_b4(8) */
      case TAIL_B4:   /* hydrocarbon_b3(7)-hydrocarbon_b4(8)-hydrocarbon_b5(9) */
	atom1 = id; 
	atom2 = id+1;
	atom3 = id+2;
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
	break;
      case ESTER_B: /* glycerol(3)-ester_B(10)-hydrocarbon_b1(11) -- SPECIAL CASE */
	atom1 = id-6; 
	atom2 = id+1;
	atom3 = id+2; 
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
	break;
      case TAIL_A2:   /* hydrocarbon_a1(5)-hydrocarbon_a2(6)-hydrocarbon_a3(7) */
      case TAIL_B2:   /* hydrocarbon_b1(5)-hydrocarbon_b2(6)-hydrocarbon_b3(7) */
	atom1 = id; 
	atom2 = id+1;
	atom3 = id+2;
	// set cisUnsat angle type for DOPC and DOPE - case for DSPC to be coded!
	fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, atom1, atom2, atom3 ); 
	break;
      default:
	printf("Error in function WriteLammpsData - default case reached. Exiting...\n"); exit(0);
      } 
    }
  }
  if (angleCount!=nAngles) {printf("Error: angleCount=%d!=%d=nAngles.\n",angleCount,nAngles); exit(0);}
}

/*******************************************************************************************************
 * DumpLammpsData -- dumps "data." file to be read into the LAMMPS program through read_data command *
 *******************************************************************************************************/
void DumpLammpsData()
{
  FILE *fp; // file pointer
  char fn[20]; // file name
  char nString[10]; 

  fn[ 0 ] = '\0';
  strcat( fn, "data." );
  
  if ( zReplicas==0 ) zReplicas = 1; // set zReplicas to default value in case it wasn't set in 'brahms.md'

  if (nLipids) {
    //    if (nDOPCsDSPCs) {
      sprintf( nString, "%d", nLipids*zReplicas );
      strcat( fn, nString );
      strcat( fn, "lip");
      //    }
  }

  if (nWaters) {
    sprintf( nString, "%d", nWaters*zReplicas );
    strcat( fn, nString );
    strcat( fn, "wat");
  }

  fp = fopen( fn, "w" );
  WriteLammpsData( fp );
  fclose(fp);
}

/**********************************************************************************************
 * WriteLammpsForcefield -- writes "forcefield" file to include into LAMMPS input script      *
 **********************************************************************************************/
void WriteLammpsForcefield( FILE *fp ) 
{
  int i,j,n;
  real springRigidity_LAMMPS_units = springRigidity * J_IN_cal / Sqr(nm_IN_Angstrom);
  real angleRigidity_LAMMPS_units = angleRigidity * J_IN_cal;
  real orientRigidity_LAMMPS_units = orientRigidity * J_IN_cal;
  real radians_in_degrees = ( 180 / PI ); // conversion factor

  //  fprintf(fp, "pair_style hybrid lj/cut 12 lj/cut/coul/cut 12 dipole/cut 9\n\n");
  fprintf(fp, "pair_style dipole/sf %2.0f\n", rCutLipLip*LENGTH_Angstrom);
  
  if (excludeBondedNebrs)
    fprintf(fp, "special_bonds lj/coul 0.0 1.0 1.0\n");
  else fprintf(fp, "special_bonds lj/coul 1.0 1.0 1.0\n");
  
  fprintf(fp, "bond_style harmonic\n");
  fprintf(fp, "angle_style hybrid cosine/squared dipole\n");
  fprintf(fp, "\n");

  for ( n = 0; n < nTypes-1; n++ ) { // write masses
    if (n==1) {
      if (nLipids==nDOPEs)
	fprintf(fp, "mass%5d%9.3f\n", n+1, massAmine);
      else if (nLipids==nDOPCsDSPCs)
	fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
      else {printf("Error: mass assignment for mixed lipid bilayers to be coded.\n"); exit(0);}
    } else
      fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
  }
  fprintf(fp, "\n");
  
  // pair_coeff for 'dipole/cut' style
  // syntax: pair_coeff iType jType ijEpsilon ijSigma cutoffLJ_optional cutoffElectro_optional
  for ( i = 0; i < nTypes-1; i++ ) { 
    for ( j = 0; j < nTypes-1; j++ ) {
      if ( i <= j ) { // LAMMPS restriction to avoid redundancy
	if ( i==1 || j==1 ){ // resolve AMINE case
	  if ( i==1 && j==1 ) {
	    if (nLipids==nDOPEs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",2,2,eps[AMINE_TYPE][AMINE_TYPE]*J_IN_cal, sig[AMINE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
	    else if (nLipids==nDOPCsDSPCs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",2,2,eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
	    else {printf("Error: pair coeff assignment for mixed lipid bilayers to be coded.\n"); exit(0);}
	  } else {
	    if (j==1) {
	      if (nLipids==nDOPEs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, 2, eps[i][AMINE_TYPE]*J_IN_cal, sig[i][AMINE_TYPE]*nm_IN_Angstrom);
	      else if (nLipids==nDOPCsDSPCs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
	      else {printf("Error: pair coeff assignment for mixed lipid bilayers to be coded.\n"); exit(0);}
	    }
	    if (i==1) {
	      if (nLipids==nDOPEs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1, eps[j][AMINE_TYPE]*J_IN_cal, sig[j][AMINE_TYPE]*nm_IN_Angstrom);
	      else if (nLipids==nDOPCsDSPCs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
	      else {printf("Error: pair coeff assignment for mixed lipid bilayers to be coded.\n"); exit(0);}
	    }
	  }
	}
	else fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
      }
    }    
  }

  fprintf(fp, "\n");

  if (nSites!=nWaters) { // to prevent segmentation fault when simulating pure water systems
    // bond_coeff for 'harmonic' style
    // syntax: bond_coeff bondType rigidityConstant referenceLength
    // NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so K_LAMMPS = k_brahms/2
    if (nLipids==nDOPCsDSPCs)
      fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", X_PHOS, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
    else if (nLipids==nDOPEs)
      fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", X_PHOS, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[AMINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
    else {printf("Error: bond assignment for mixed lipid bilayers to be coded.\n"); exit(0);}
    fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", PHOS_GLY, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", GLY_EST, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", EST_CHAIN, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d  %7.3f%7.3f\n", CHAIN_CHAIN, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
    //  fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14, springRigidity_LAMMPS_units/2, scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
    fprintf(fp, "\n");

    // angle_coeff for 'cosine/squared' style
    // syntax: angle_coeff angleType rigidityConstant referenceAngle
    // NOTE: in LAMMPS, this potential is defined as K[cosTheta-cosTheta0]^2, so K_LAMMPS = k_brahms/2
    fprintf(fp, "angle_coeff%5d  cosine/squared%8.3f%8.2f\n", _phos_, angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d  cosine/squared%8.3f%8.2f\n", _gly_, angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d  cosine/squared%8.3f%8.2f\n", _piSat_, angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d  cosine/squared%8.3f%8.2f\n", _cisUnsat_, angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d  dipole        %8.3f%8.2f\n", DIPOLE, orientRigidity_LAMMPS_units/2, refAngle_orient*radians_in_degrees); 
    fprintf(fp, "\n");
  }
}

/*******************************************************************************************
 * DumpLammpsForcefield -- dumps "forcefield." file to be included into the LAMMPS program *
 *******************************************************************************************/
void DumpLammpsForcefield()
{
  FILE *fp; // file pointer
  char fn[20]; // file name

  fn[ 0 ] = '\0';
  strcat( fn, "forcefield.elba" );
  
  fp = fopen( fn, "w" );
  WriteLammpsForcefield( fp );
  fclose(fp);
}
